/*
 * Copyright 2011-2013 Blender Foundation
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "stdosl.h"
#include "node_texture.h"

/* Musgrave fBm
 *
 * H: fractal increment parameter
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 *
 * from "Texturing and Modelling: A procedural approach"
 */

float noise_musgrave_fBm(point ip, float H, float lacunarity, float octaves)
{
  float rmd;
  float value = 0.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);
  int i;
  point p = ip;

  for (i = 0; i < (int)octaves; i++) {
    value += safe_noise(p, "signed") * pwr;
    pwr *= pwHL;
    p *= lacunarity;
  }

  rmd = octaves - floor(octaves);
  if (rmd != 0.0)
    value += rmd * safe_noise(p, "signed") * pwr;

  return value;
}

/* Musgrave Multifractal
 *
 * H: highest fractal dimension
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 */

float noise_musgrave_multi_fractal(point ip, float H, float lacunarity, float octaves)
{
  float rmd;
  float value = 1.0;
  float pwr = 1.0;
  float pwHL = pow(lacunarity, -H);
  int i;
  point p = ip;

  for (i = 0; i < (int)octaves; i++) {
    value *= (pwr * safe_noise(p, "signed") + 1.0);
    pwr *= pwHL;
    p *= lacunarity;
  }

  rmd = octaves - floor(octaves);
  if (rmd != 0.0)
    value *= (rmd * pwr * safe_noise(p, "signed") + 1.0); /* correct? */

  return value;
}

/* Musgrave Heterogeneous Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hetero_terrain(
    point ip, float H, float lacunarity, float octaves, float offset)
{
  float value, increment, rmd;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;
  int i;
  point p = ip;

  /* first unscaled octave of function; later octaves are scaled */
  value = offset + safe_noise(p, "signed");
  p *= lacunarity;

  for (i = 1; i < (int)octaves; i++) {
    increment = (safe_noise(p, "signed") + offset) * pwr * value;
    value += increment;
    pwr *= pwHL;
    p *= lacunarity;
  }

  rmd = octaves - floor(octaves);
  if (rmd != 0.0) {
    increment = (safe_noise(p, "signed") + offset) * pwr * value;
    value += rmd * increment;
  }

  return value;
}

/* Hybrid Additive/Multiplicative Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_hybrid_multi_fractal(
    point ip, float H, float lacunarity, float octaves, float offset, float gain)
{
  float result, signal, weight, rmd;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;
  int i;
  point p = ip;

  result = safe_noise(p, "signed") + offset;
  weight = gain * result;
  p *= lacunarity;

  for (i = 1; (weight > 0.001) && (i < (int)octaves); i++) {
    if (weight > 1.0)
      weight = 1.0;

    signal = (safe_noise(p, "signed") + offset) * pwr;
    pwr *= pwHL;
    result += weight * signal;
    weight *= gain * signal;
    p *= lacunarity;
  }

  rmd = octaves - floor(octaves);
  if (rmd != 0.0)
    result += rmd * ((safe_noise(p, "signed") + offset) * pwr);

  return result;
}

/* Ridged Multifractal Terrain
 *
 * H: fractal dimension of the roughest area
 * lacunarity: gap between successive frequencies
 * octaves: number of frequencies in the fBm
 * offset: raises the terrain from `sea level'
 */

float noise_musgrave_ridged_multi_fractal(
    point ip, float H, float lacunarity, float octaves, float offset, float gain)
{
  float result, signal, weight;
  float pwHL = pow(lacunarity, -H);
  float pwr = pwHL;
  int i;
  point p = ip;

  signal = offset - fabs(safe_noise(p, "signed"));
  signal *= signal;
  result = signal;
  weight = 1.0;

  for (i = 1; i < (int)octaves; i++) {
    p *= lacunarity;
    weight = clamp(signal * gain, 0.0, 1.0);
    signal = offset - fabs(safe_noise(p, "signed"));
    signal *= signal;
    signal *= weight;
    result += signal * pwr;
    pwr *= pwHL;
  }

  return result;
}

/* Shader */

shader node_musgrave_texture(
    int use_mapping = 0,
    matrix mapping = matrix(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    string type = "fBM",
    float Dimension = 2.0,
    float Lacunarity = 1.0,
    float Detail = 2.0,
    float Offset = 0.0,
    float Gain = 1.0,
    float Scale = 5.0,
    point Vector = P,
    output float Fac = 0.0,
    output color Color = 0.0)
{
  float dimension = max(Dimension, 1e-5);
  float octaves = clamp(Detail, 0.0, 16.0);
  float lacunarity = max(Lacunarity, 1e-5);
  float intensity = 1.0;

  point p = Vector;

  if (use_mapping)
    p = transform(mapping, p);

  p = p * Scale;

  if (type == "multifractal")
    Fac = intensity * noise_musgrave_multi_fractal(p, dimension, lacunarity, octaves);
  else if (type == "fBM")
    Fac = intensity * noise_musgrave_fBm(p, dimension, lacunarity, octaves);
  else if (type == "hybrid_multifractal")
    Fac = intensity *
          noise_musgrave_hybrid_multi_fractal(p, dimension, lacunarity, octaves, Offset, Gain);
  else if (type == "ridged_multifractal")
    Fac = intensity *
          noise_musgrave_ridged_multi_fractal(p, dimension, lacunarity, octaves, Offset, Gain);
  else if (type == "hetero_terrain")
    Fac = intensity * noise_musgrave_hetero_terrain(p, dimension, lacunarity, octaves, Offset);

  Color = color(Fac, Fac, Fac);
}
